FUNCTION delta(a, r)
    IMPLICIT NONE
    DOUBLE PRECISION Delta, a, r
    delta = r**2 + a**2 - 2. * r
END

FUNCTION P(a, lam, r)
    IMPLICIT NONE
    DOUBLE PRECISION P, a, lam, r
    P = r**2 + a**2 - lam * a
END

! accordingt to Fanton 1997, to make sure the photon can escape
FUNCTION RL(a, lam, r)
    IMPLICIT NONE
    DOUBLE PRECISION RL, a, lam, r, delta
    RL = ((2d0 * r - a * lam)**2 + delta(a, r) * (r**2 + 2d0 * r)) / delta(a, r)
END


FUNCTION RTHETA(a, lam, q, u)
    IMPLICIT NONE
    DOUBLE PRECISION RTHETA, a, lam, q, theta, u
    RTHETA = q**2 + (a * u)**2 - (lam * u)**2 / (1d0 - u**2)
    RETURN
END

FUNCTION RR(a, lam, q, r)
    IMPLICIT NONE
    DOUBLE PRECISION RR, a, lam, q, r, delta, P
    RR = P(a, lam, r)**2 - delta(a, r) * (q**2 + (lam - a)**2)
    RETURN
END


FUNCTION TrajR(a, lam, q, r)
    IMPLICIT NONE
    DOUBLE PRECISION TrajR, a, lam, q, r, RR, P, delta
    TrajR = 1d0 / Sqrt(RR(a, lam, q, r))
    RETURN
END

! integration over theta need to change to du, which leads to -sin(theta)
! Fortunately, theta is positive and smaller than pi/2, which means there is no need to consider the sign.
FUNCTION TrajTheta(a, lam, q, u)
    IMPLICIT NONE
    DOUBLE PRECISION TrajTheta, a, lam, q, theta, RTHETA, u
    TrajTheta = 1d0 / Sqrt(RTHETA(a, lam, q, u) * (1d0 - u**2))
    RETURN
END


FUNCTION upositive(a, lam, q)
    IMPLICIT NONE
    DOUBLE PRECISION upositive, a, lam, q, xt
    xt = lam**2 + q**2 - a**2
    upositive = (sqrt(xt**2 + 4d0 * (a * q)**2) - xt) / (2d0 * a)
    RETURN
END


FUNCTION phir(a, lam, q, r)
    IMPLICIT NONE
    DOUBLE PRECISION phir, a, lam, q, r, RR, P, delta
    phir = (2.d0 * a * r - lam * a**2) / Sqrt(RR(a, lam, q, r)) / delta(a, r)
    RETURN
END

FUNCTION phitheta(a, lam, q, u)
    IMPLICIT NONE
    DOUBLE PRECISION phitheta, a, lam, q, theta, RTHETA, u
    phitheta = lam / (1d0 - u**2) / Sqrt(RTHETA(a, lam, q, u) * (1d0 - u**2))
    RETURN
END

FUNCTION tr(a, lam, q, r)
    IMPLICIT NONE
    DOUBLE PRECISION tr, a, lam, q, r, RR, P, delta
    tr = (r**2 * (r**2 + a**2) + 2.d0 * a * r * (a - lam)) / Sqrt(RR(a, lam, q, r))
    tr = tr / delta(a, r)
    RETURN
END

FUNCTION ttheta(a, lam, q, u)
    IMPLICIT NONE
    DOUBLE PRECISION ttheta, a, lam, q, theta, RTHETA, u
    ttheta = a**2 * u**2 / Sqrt(RTHETA(a, lam, q, u) * (1d0 - u**2))
    RETURN
END


PROGRAM main
    IMPLICIT NONE
    DOUBLE PRECISION delta, P, RTHETA, RR, phir, xr, xth, &
            phitheta, tr, ttheta, sr, pi, stheta, xr1, xr2, xth1, xth2&
            , RTHETA2, RR2, phir2, phitheta2, tr2, ttheta2, RL, xrl&
            , TrajR, TrajTheta, strar, stratheta, u, up, u1, u2&
            , sphiu1, sphiu2, stu1, stu2, upositive, strth
    DOUBLE PRECISION a0, lam0, q0, th1, th2, r1, r2, sphi, eps
    DOUBLE PRECISION rx, thetax, x, y, r, g, ds, str, stth
    INTEGER    i, j, k, maxn, rin, rout, a_step, th2_step
    CHARACTER (len = 30)  filenamein, filenameout
    CHARACTER (len = 4) char_a0,char_th2,char_r2
    EXTERNAL phir, phitheta, tr, ttheta, TrajR, TrajTheta

    Real time_begin, time_end

    pi = 4d0 * Atan(1d0)
    th1 = pi / 2.d0
    u1 = 0d0
    eps = 1d-2
    maxn = 80000000



    !设置参数
    !获取命令参数
    CALL GET_COMMAND_ARGUMENT(1 , char_a0)
    CALL GET_COMMAND_ARGUMENT(2 , char_th2)
    CALL GET_COMMAND_ARGUMENT(3 , char_r2)
    CALL GET_COMMAND_ARGUMENT(4 , filenamein)
    CALL GET_COMMAND_ARGUMENT(5 , filenameout)

    read(char_a0,*) a0
    read(char_th2,*) th2
    read(char_r2,*) r2


    th2 = th2 * pi / 180d0
    u2 = cos(th2) * 0.9999d0
    OPEN(22, file = filenamein, status = 'old')
    OPEN(33, file = filenameout, status = 'unknown')

    DO i = 1, maxn
        read(22, *, end = 1026) x, y, r1, g, ds, lam0, q0
        q0 = sqrt(q0)
        xrl = RL(a0, lam0, r1) - lam0**2 - q0**2
        IF(xrl < 0d0) THEN
            WRITE(333, *) "L>LR", x, y, r1, g, ds, lam0, q0
            GOTO 666
        ENDIF

        xr = RR(a0, lam0, q0, r1)
        xth = RTHETA(a0, lam0, q0, u2)
        up = upositive(a0, lam0, q0)


        IF(xr >= 0d0 .and. xth >= 0d0) THEN
            CALL simp(TrajR, a0, lam0, q0, r1, r2, strar)
            CALL qromb(TrajTheta, a0, lam0, q0, u1, up, sphiu1)
            CALL qromb(TrajTheta, a0, lam0, q0, u2, up, sphiu2)
        ELSE
            CALL simp(TrajR, a0, lam0, q0, r1 + eps, r2, strar)
            CALL qromb(TrajTheta, a0, lam0, q0, u1 + 0.0001d0 * eps, up, sphiu1)
            CALL qromb(TrajTheta, a0, lam0, q0, u2 + 0.0001d0 * eps, up, sphiu2)

        ENDIF

        IF(y <0d0) THEN
            strth = sphiu1 - sphiu2
        ELSE
            strth = sphiu1 + sphiu2
        ENDIF

        strar = strar - strth

        IF(xr >= 0d0) THEN
            CALL simp(phir, a0, lam0, q0, r1, r2, sr)
            CALL simp(tr, a0, lam0, q0, r1, r2, str)
        ELSE
            !to remove the cases of negative R due to precision of calculation from CXX file.
            CALL simp(phir, a0, lam0, q0, r1 + eps, r2, sr)
            CALL simp(tr, a0, lam0, q0, r1 + eps, r2, str)
        ENDIF

        IF(xth >= 0d0) THEN
            CALL qromb(phitheta, a0, lam0, q0, u1, up, sphiu1)
            CALL qromb(ttheta, a0, lam0, q0, u1, up, stu1)
            CALL qromb(phitheta, a0, lam0, q0, u2, up, sphiu2)
            CALL qromb(ttheta, a0, lam0, q0, u2, up, stu2)
        ELSE
            CALL qromb(phitheta, a0, lam0, q0, u1, up, sphiu1)
            CALL qromb(ttheta, a0, lam0, q0, u1, up, stu1)
            CALL qromb(phitheta, a0, lam0, q0, u2 + 0.01d0 * eps, up, sphiu2)
            CALL qromb(ttheta, a0, lam0, q0, u2 + 0.01d0 * eps, up, stu2)
        ENDIF

        IF(y <0d0) THEN
            stheta = sphiu1 - sphiu2
            stth = stu1 - stu2
        ELSE
            stheta = sphiu1 + sphiu2
            stth = stu1 + stu2
        ENDIF
        write(33, *) x, y, r1, g, ds, lam0, q0, sr, stheta, str, stth, strar, strth
        666    CONTINUE
    ENDDO

    1026    CONTINUE
    CLOSE(22)
    CLOSE(33)

END


SUBROUTINE qromb(func, spin, slam, sq, a, b, ss)
    INTEGER JMAX, JMAXP, K, KM
    DOUBLE PRECISION a, b, func, ss, EPS, spin, slam, sq
    EXTERNAL func
    PARAMETER (EPS = 1.d-6, JMAX = 25, JMAXP = JMAX + 1, K = 5, KM = K - 1)
    !    CU    USES polint,trapzd
    INTEGER j
    DOUBLE PRECISION dss, h(JMAXP), s(JMAXP)
    h(1) = 1.d0
    do j = 1, JMAX
        call trapzd(func, spin, slam, sq, a, b, s(j), j)
        if (j>=K) then
            call polint(h(j - KM), s(j - KM), K, 0.d0, ss, dss)
            if (abs(dss)<=EPS * abs(ss)) then
                return
            end if
        endif
        s(j + 1) = s(j)
        h(j + 1) = 0.25d0 * h(j)
    end do
    write(*, *) 'too many steps in qromb', spin, slam, sq, a, b
    ss = 0.0
END
!C  (C) Copr. 1986-92 Numerical Recipes Software &H1216.



SUBROUTINE trapzd(func, spin, slam, sq, a, b, s, n)
    INTEGER n
    DOUBLE PRECISION a, b, s, func, spin, slam, sq
    EXTERNAL func
    INTEGER it, j
    DOUBLE PRECISION del, sum, tnm, x
    if (n==1) then
        s = 0.5d0 * (b - a) * (func(spin, slam, sq, a) + func(spin, slam, sq, b))
    else
        it = 2**(n - 2)
        tnm = it
        del = (b - a) / tnm
        x = a + 0.5d0 * del
        sum = 0.d0
        do j = 1, it
            sum = sum + func(spin, slam, sq, x)
            x = x + del
        end do
        s = 0.5d0 * (s + (b - a) * sum / tnm)
    endif
    return
END
!C  (C) Copr. 1986-92 Numerical Recipes Software &H1216.


SUBROUTINE polint(xa, ya, n, x, y, dy)
    INTEGER n, NMAX
    DOUBLE PRECISION dy, x, y, xa(n), ya(n)
    PARAMETER (NMAX = 10)
    INTEGER i, m, ns
    DOUBLE PRECISION den, dif, dift, ho, hp, w, c(NMAX), d(NMAX)
    ns = 1
    dif = abs(x - xa(1))
    do i = 1, n
        dift = abs(x - xa(i))
        if (dift<dif) then
            ns = i
            dif = dift
        endif
        c(i) = ya(i)
        d(i) = ya(i)
    end do
    y = ya(ns)
    ns = ns - 1
    do m = 1, n - 1
        do i = 1, n - m
            ho = xa(i) - x
            hp = xa(i + m) - x
            w = c(i + 1) - d(i)
            den = ho - hp
            if(den==0.d0)write(*, *) 'failure in polint'
            den = w / den
            d(i) = hp * den
            c(i) = ho * den
        end do
        if (2 * ns<n - m)then
            dy = c(ns + 1)
        else
            dy = d(ns)
            ns = ns - 1
        endif
        y = y + dy
    end do
    return
END
!C  (C) Copr. 1986-92 Numerical Recipes Software &H1216.

subroutine simp(func, spin, slam, sq, a, b, ss)
    IMPLICIT NONE
    DOUBLE PRECISION func, spin, slam, sq, a, b, ss, epslion, s0
    EXTERNAL func
    epslion = 1d-6
    s0 = 0d0
    call adaptiveSim(func, spin, slam, sq, a, b, epslion, s0, ss)
end

recursive SUBROUTINE adaptiveSim(func, spin, slam, sq, a, b, epsilon, S, res)
    IMPLICIT NONE
    DOUBLE PRECISION  func, spin, slam, sq, a, b, epsilon, S, res, h, S1, S2, res1, res2, RTHETA
    EXTERNAL func, RTHETA
    h = 5d-1 * (b - a)
    S1 = h / 6 * (func(spin, slam, sq, a) + 4 * func(spin, slam, sq, a + h / 2) + func(spin, slam, sq, a + h))
    S2 = h / 6 * (func(spin, slam, sq, a + h) + 4 * func(spin, slam, sq, a + 3 * h / 2) + func(spin, slam, sq, b))

    if (abs(S - S1 - S2) <= epsilon) then
        res = (S1 + S2) + (S1 + S2 - S) / 15
    else
        call adaptiveSim(func, spin, slam, sq, a, (a + b) / 2, epsilon, S1, res1)
        call adaptiveSim(func, spin, slam, sq, (a + b) / 2, b, epsilon, S2, res2)
        res = res1 + res2
    end if
    RETURN
end